In 2017, the Tuveson Lab at Cold Spring Harbor Cancer Center published a paper written by Elyada et al that detailed the discovery of cancer-associated fibroblasts (CAFs) in mice. The subtypes were then validated in human samples affected with PDAC in a subsequent paper released in 2019. Here we will use SCISSORS to identify the CAF subtypes within the larger stroma population, fine-grained immune cell types within broadly-defined immune clusters, and ductal & PDAC subtypes within the ductal group. You can install SCISSORS from our GitHub repository.
library(VAM) # single cell GSEA
library(dplyr) # tidy data
library(Seurat) # single cell infrastructure
library(ggplot2) # pretty plots
library(SingleR) # cell type assignment
library(decoderr) # de novo deconvolution
library(SCISSORS) # our package
library(mixtools) # Gaussian mixture model estimation
library(patchwork) # align plots
library(latex2exp) # LaTeX
library(paletteer) # color palettes
library(CONICSmat) # CNV estimation
library(reticulate) # Python interface
library(kableExtra) # pretty tables
library(wesanderson) # more color palettesimport numpy as np
from openTSNE import TSNEEmbedding
from openTSNE import initialization
from openTSNE.affinity import Multiscale
from openTSNE.affinity import PerplexityBasedNNFirst we load in the \(\text{gene} \times \text{cell}\) counts matrix, then create a Seurat object to hold it in. Next, we add samplename, tissue type, and patient sex metadata taken from the publicly available dataset.
raw_counts <- Read10X(data.dir = "~/Desktop/Data/Elyada Raw/All Human/")
pdac <- CreateSeuratObject(raw_counts,
project = "Elyada",
min.cells = 3,
min.features = 500)
pdac@meta.data$sample <- case_when(grepl("-1", rownames(pdac@meta.data)) ~ "SRR9274536",
grepl("-2", rownames(pdac@meta.data)) ~ "SRR9274537",
grepl("-3", rownames(pdac@meta.data)) ~ "SRR9274538",
grepl("-4", rownames(pdac@meta.data)) ~ "SRR9274539",
grepl("-5", rownames(pdac@meta.data)) ~ "SRR9274540",
grepl("-6", rownames(pdac@meta.data)) ~ "SRR9274541",
grepl("-7", rownames(pdac@meta.data)) ~ "SRR9274542",
grepl("-8", rownames(pdac@meta.data)) ~ "SRR9274543",
grepl("-9", rownames(pdac@meta.data)) ~ "SRR9274544")
pdac@meta.data$condition <- case_when(grepl("-1", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-2", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-3", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-4", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-5", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-6", rownames(pdac@meta.data)) ~ "AdjNorm",
grepl("-7", rownames(pdac@meta.data)) ~ "PDAC",
grepl("-8", rownames(pdac@meta.data)) ~ "AdjNorm",
grepl("-9", rownames(pdac@meta.data)) ~ "PDAC")
pdac@meta.data$sex <- case_when(grepl("-1", rownames(pdac@meta.data)) ~ "female",
grepl("-2", rownames(pdac@meta.data)) ~ "male",
grepl("-3", rownames(pdac@meta.data)) ~ "male",
grepl("-4", rownames(pdac@meta.data)) ~ "male",
grepl("-5", rownames(pdac@meta.data)) ~ "male",
grepl("-6", rownames(pdac@meta.data)) ~ "female",
grepl("-7", rownames(pdac@meta.data)) ~ "female",
grepl("-8", rownames(pdac@meta.data)) ~ "male",
grepl("-9", rownames(pdac@meta.data)) ~ "female")We run the typical single cell pre-processing steps on our cells - normalization, dimension reduction, and clustering.
pdac <- PrepareData(seurat.object = pdac,
n.HVG = 4000,
n.PC = 20,
regress.mt = FALSE,
regress.cc = FALSE,
which.dim.reduc = "tsne",
initial.resolution = 0.4,
random.seed = 629)## [1] "Running t-SNE on 20 principal components with perplexity = 30"
## [1] "Found 12 unique clusters"
Let’s check out the t-SNE embedding. It looks decent, but could definitely be improved. Globally, the clusters are arranged in a blob - which isn’t very informative - though the local structure seems to have been preserved fairly well. Visually, there are several clusters that contain subgroups - clusters 3, 4, 5, 6, & 9 look like good candidates for reclustering.
p0a <- DimPlot(pdac, reduction = "tsne", pt.size = 0.75) +
scale_color_paletteer_d("ggthemes::Classic_20") +
labs(x = "t-SNE 1", y = "t-SNE 2") +
theme_yehlab() +
theme(legend.text = element_text(size = 10)) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p0aWe compute the distribution of the silhouette scores for each cluster.
sil_df <- ComputeSilhouetteScores(pdac, avg = FALSE)
p0b <- ggplot(sil_df, aes(x = Cluster, y = Score, fill = Cluster)) +
geom_violin(draw_quantiles = .5, color = "black", scale = "width") +
scale_fill_paletteer_d("ggthemes::Classic_20") +
labs(y = "Silhouette Score", fill = "Louvain\nClusters") +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(fill = NA, size = 1),
legend.position = "right",
legend.direction = "vertical",
axis.title.x = element_blank(),
legend.title = element_text(size = 22),
legend.text = element_text(size = 16),
axis.title.y = element_text(size = 22),
axis.text = element_text(size = 16),
axis.ticks = element_line(),
plot.subtitle = element_text(face = "italic", size = 10))Looking at the silhouette score distribution for each cluster, we see that Cluster 10 seems to have the best fit, and Clusters 1, 2, & 9 seem to have pretty poor fits.
p0bLastly we’ll identify marker genes for each of the clusters.
pdac_markers <- FindAllMarkers(pdac,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p0c <- DotPlot(pdac, features = unique(pdac_markers$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed", y = "Louvain Cluster") +
theme(axis.text.x = element_text(angle = 90, size = 10, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 14),
axis.text.y = element_text(size = 18)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p0cI think the two-dimensional visualization of the cells could be improved. We’ll try using UMAP and the Fast Fourier Transform-accelerated Fit-SNE (as implemented in the openTSNE library) to improve the embedding.
It seems like UMAP does a good job of clearly separating our clusters and preserving the global structure of the data. However, it’s difficult to see local structure within some of the clusters due to their density.
pdac <- RunUMAP(pdac,
reduction = "pca",
dims = 1:20,
umap.method = "uwot",
n.components = 2,
n.epochs = 750,
n.neighbors = 50,
metric = "cosine",
seed.use = 629,
verbose = FALSE)
p1 <- DimPlot(pdac, reduction = "umap", pt.size = 0.75) +
scale_color_paletteer_d("ggthemes::Classic_20") +
labs(x = "UMAP 1", y = "UMAP 2") +
theme_yehlab() +
theme(legend.text = element_text(size = 10)) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p1For information on how to install the openTSNE implementation of Fit-SNE and how to run the algorithm, please visit the excellent GitHub repository of Pavlin Policar.
First we’ll run a simple, standard Fit-SNE embedding. It’s necessary to make the PCA embeddings accessible by Python.
pc_df <- Embeddings(pdac, reduction = "pca")# import data
pc_df = np.array(r.pc_df)
# run Fit-SNE
affin = PerplexityBasedNN(pc_df, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(pc_df, random_state=629)
tsne1 = TSNEEmbedding(init, affin, negative_gradient_method='fft')
embed1 = tsne1.optimize(n_iter=350, exaggeration=12, momentum=0.6)
embed2 = embed1.optimize(n_iter=750, momentum=0.8)The embedding looks good, so we’ll use it going forwards.
embed <- as.matrix(py$embed2)
rownames(embed) <- colnames(pdac)
pdac@reductions$fitsne <- CreateDimReducObject(embeddings = embed,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p2 <- DimPlot(pdac, reduction = "fitsne", pt.size = 0.75) +
scale_color_paletteer_d("ggthemes::Classic_20") +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(legend.text = element_text(size = 10)) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p2Due to the way Seurat accesses cell embeddings, we’ll need to replace our original t-SNE dimension reduction in our Seurat object with the new Fit-SNE version. We’ll keep the original Barnes-Hut t-SNE embedding under a separate name.
pdac@reductions$bh_tsne <- pdac@reductions$tsne
pdac@reductions$tsne <- pdac@reductions$fitsneHere we use the SingleR package to identify broad cell types. The reference dataset we load is an sctransform-normalized version of the raw counts available in scRNAseq::BaronPancreasData(), which consists of normal pancreas cells that were sequenced and annotated by the researchers.
sc_ref <- readRDS("/Volumes/labs/Home/Jen Jen Yeh Lab/Jack/scRNAseq/Seurat/single_cell_ref_normalized.Rds")
sc_preds <- SingleR(test = data.frame(pdac@assays$SCT@data),
ref = sc_ref,
labels = sc_ref$label,
method = "cluster",
clusters = pdac$seurat_clusters,
de.method = "wilcox")
pdac[["SingleR.labels.sc"]] <- sc_preds$labels[match(pdac[[]][["seurat_clusters"]], rownames(sc_preds))]
pdac$SingleR.labels.sc <- case_when(pdac$SingleR.labels.sc == "acinar" ~ "Acinar",
pdac$SingleR.labels.sc == "activated_stellate" ~ "Activated Stellate",
pdac$SingleR.labels.sc == "ductal" ~ "Ductal",
pdac$SingleR.labels.sc == "macrophage" ~ "Macrophage",
pdac$SingleR.labels.sc == "t_cell" ~ "T")We can see that there’s a large immune population, as well as smaller ductal, fibroblast (denoted activated stellate in the reference dataset), and acinar groups. These broad cell types line up with what we expected to see given the authors’ original cell cluster annotations.
p3 <- DimPlot(pdac, reduction = "tsne", group.by = "SingleR.labels.sc", pt.size = 0.75) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 4)))
p3This dataset is composed of labeled & log-normalized bulk RNA-seq samples from the Human Primary Cell Atlas.
bulk_ref <- HumanPrimaryCellAtlasData()
bulk_preds <- SingleR(test = data.frame(pdac@assays$SCT@data),
ref = bulk_ref,
labels = bulk_ref$label.main,
method = "cluster",
clusters = pdac$seurat_clusters,
de.method = "wilcox")
pdac[["SingleR.labels.bulk"]] <- bulk_preds$labels[match(pdac[[]][["seurat_clusters"]],
rownames(bulk_preds))]
pdac$SingleR.labels.bulk <- case_when(pdac$SingleR.labels.bulk == "B_cell" ~ "B",
pdac$SingleR.labels.bulk == "Epithelial_cells" ~ "Epithelial",
pdac$SingleR.labels.bulk == "B_cell-" ~ "B",
pdac$SingleR.labels.bulk == "T_cells" ~ "T",
TRUE ~ pdac$SingleR.labels.bulk)The bulk reference gives us somewhat more granular labels for the immune cells, and confirms the identities of the ductal / epithelial and stroma clusters.
p4 <- DimPlot(pdac, reduction = "tsne", group.by = "SingleR.labels.bulk", pt.size = 0.75) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 4)))
p4One shouldn’t use SingleR as the final authority for cell types, but we were able to confirm the identities of the ductal and fibroblast clusters, which is important for this dataset as the cells we are most interested in are the cancer-associated fibroblasts (CAFs).
Next we’ll attempt to identify malignant cells using single-cell copy number variation estimation as implemented in the CONCISmat package. Details of the GMM methodology used can be found at the Diaz Lab’s GitHub repository. Note: this step is memory-intensive because 1) it requires the sparse counts matrix to be cast to a dense matrix and 2) a lot of Gaussian mixture models get estimated. If your machine doesn’t have a lot of RAM it might be best to skip this and manually annotate the malignant cells instead through the usage of canonical markers or the VAM single-cell GSEA methodology.
chrom_regions <- read.table("/Volumes/labs/Home/Jen Jen Yeh Lab/Jack/scRNAseq/chrom_arm_positions.txt",
sep = "\t",
row.names = 1,
header = TRUE)
gene_pos <- getGenePositions(rownames(pdac))
cpm <- t(t(as.matrix(pdac@assays$SCT@counts)) / colSums(as.matrix(pdac@assays$SCT@counts))) * 10^5
cpm <- log2(cpm + 1)
norm_factor <- calcNormFactors(cpm)
cnv_est <- plotAll(mat = cpm,
normFactor = norm_factor,
regions = chrom_regions,
gene_pos = gene_pos,
fname = "Elyada")After estimating CNVs, we cluster the cells into \(k = 3\) clusters, with the hope of finding one large cluster of normal cells and two smaller clusters composed of CAFs and PDAC cells.
bic_table <- read.table("./Elyada_BIC_LR.txt",
sep = "\t",
row.names = 1,
header = TRUE,
check.names = FALSE)
cand_regions <- rownames(bic_table[bic_table$`BIC difference` > 1000 & bic_table$`LRT adj. p-val` < .01, ])
hist1 <- plotHistogram(cnv_est[, cand_regions],
cpm,
clusters = 3,
zscoreThreshold = 3,
celltypes = pdac$SingleR.labels.bulk,
patients = pdac$sample)We add the normal vs. malignant cell labels in to our Seurat object’s metadata, then visualize the results. As expected, the malignant cells are located in the clusters identified by SingleR as ductal cells and fibroblasts. This indicates that CONICSmat did a solid job of estimating the CNVs - no easy feat with sparse, noisy single-cell data.
normal <- which(hist1 == 1)
malignant <- which(hist1 != 1)
pdac@meta.data$malig <- ifelse(rownames(pdac@meta.data) %in% names(normal), "Normal", "Malignant")
p5 <- DimPlot(pdac, reduction = "tsne", group.by = "malig", pt.size = 0.75) +
scale_color_manual(values = wes_palette("Zissou1", n = 5)[c(5, 2)]) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p5Next, we use Dr. Xianlu Peng’s DECODER in order to deconvolve the dataset and assign weights to each cell using non-negative matrix factorization. We calculate basal & classical PDAC, normal & activated stroma, immune, and endocrine & exocrine pancreas compartment weights.
sample_wts_unscaled <- Decon_single_sample("TCGA_RNAseq_PAAD", pdac@assays$SCT@data, "geneSymbol")
sample_wts <- Norm_PDAC_weights(sample_wts_unscaled)
pdac <- AddMetaData(pdac, sample_wts$Immune, "immune")
pdac <- AddMetaData(pdac, sample_wts$bcRatio, "bc_ratio")
pdac <- AddMetaData(pdac, sample_wts$Exocrine, "exocrine")
pdac <- AddMetaData(pdac, sample_wts$Endocrine, "endocrine")
pdac <- AddMetaData(pdac, sample_wts_unscaled[, 9], "basal")
pdac <- AddMetaData(pdac, sample_wts_unscaled[, 5], "classical")
pdac <- AddMetaData(pdac, sample_wts$NormalStroma, "norm_stroma")
pdac <- AddMetaData(pdac, sample_wts$ActivatedStroma, "act_stroma")The basal weights are highest in a subcluster of the ductal group identified through SingleR. This is interesting as the authors did not find evidence of the basal subtype in their paper.
p6 <- FeaturePlot(pdac, reduction = "tsne", features = "basal", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p6The classical weights are highest in another subcluster of the ductal cluster, and cells with high classical weights do not collocate with those that have high basal weights. The putative classical and basal PDAC cells also align closely with the cells identified through CONCISmat as malignant.
p7 <- FeaturePlot(pdac, reduction = "tsne", features = "classical", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p7The cluster identified through SingleR as acinar cells is the only cluster with high exocrine pancreas weights.
p8 <- FeaturePlot(pdac, reduction = "tsne", features = "exocrine", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p8No cells have high endocrine pancreas weights.
p9 <- FeaturePlot(pdac, reduction = "tsne", features = "endocrine", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p9Once again, we confirm the largeness of the immune cell population in this dataset.
p10 <- FeaturePlot(pdac, reduction = "tsne", features = "immune", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p10Cells with high normal stroma stroma weights are located in the cluster identified by SingleR as being stromal cells.
p11 <- FeaturePlot(pdac, reduction = "tsne", features = "norm_stroma", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p11Cells with high activated stroma weights are also located in the fibroblast cluster, and do not intersect with the cells that have high normal stroma scores. This indicates that SCISSORS will most likely perform well on the fibroblast cluster and be able to quickly tease out the cell subtypes.
p12 <- FeaturePlot(pdac, reduction = "tsne", features = "act_stroma", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme(plot.title = element_blank()) +
theme_yehlab() +
NoLegend()
p12Now that we have rough labels from SingleR, CNVs from CONICSmat, and compartment weights from DECODER, we should have more than enough cell-level metadata to look for and annotate cell subtypes using SCISSORS.
The fibroblast marker genes provided by Elyada et al match the SingleR results defining cluster 6 as containing fibroblasts.
p13 <- FeaturePlot(pdac, reduction = "tsne", features = "COL1A1", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "COL1A1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p14 <- FeaturePlot(pdac, reduction = "tsne", features = "COL3A1", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "COL3A1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p15 <- FeaturePlot(pdac, reduction = "tsne", features = "LUM", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "LUM") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p16 <- FeaturePlot(pdac, reduction = "tsne", features = "DCN", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "DCN") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p13 | p14) / (p15 | p16)Here we use ReclusterCells() to identify subclusters within cluster 6. We find five distinct subclusters.
fibro <- ReclusterCells(pdac,
which.clust = 6,
n.HVG = 4000,
n.PC = 20,
resolution.vals = c(.03, .05, .1),
k.vals = c(20, 30, 40),
redo.embedding = TRUE)
fibro_pc <- Embeddings(fibro, "pca")## [1] "Reclustering cells in cluster 6 using k = 20 & resolution = 0.03; S = 0.579"
We’ll again run Fit-SNE on the reclustered cells, for consistencies sake.
# import data
fibro_pc = r.fibro_pc
# run Fit-SNE
affin_fibro = PerplexityBasedNN(fibro_pc, perplexity=40, metric='cosine', random_state=629)
init = initialization.pca(fibro_pc, random_state=629)
tsne_f = TSNEEmbedding(init, affin_fibro, negative_gradient_method='fft')
embed_f1 = tsne_f.optimize(n_iter=250, exaggeration=5, momentum=0.4)
embed_f2 = embed_f1.optimize(n_iter=750, exaggeration=1, momentum=0.8)Pulling the results back into R and visualizing them shows clear separation between the subclusters. There’s some noise in subcluster 0, but other than that the reembedding looks solid.
embed_fibro <- as.matrix(py$embed_f2)
rownames(embed_fibro) <- colnames(fibro)
fibro@reductions$bh_tsne <- fibro@reductions$tsne
fibro@reductions$tsne<- CreateDimReducObject(embeddings = embed_fibro,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p17 <- DimPlot(fibro, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p17First we identify the endothelial and perivascular cells using PLVAP and RGS5 expression.
p18 <- FeaturePlot(fibro, reduction = "tsne", features = "PLVAP", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "PLVAP") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p19 <- FeaturePlot(fibro, reduction = "tsne", features = "RGS5", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "RGS5") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p18 | p19) / p17Next, we use the VAM method of single cell gene set enrichment analysis to determine which clusters are enriched for the iCAF and myCAF marker genes, as well as the general pan-CAF marker set. We use the marker genes identified by the authors.
icaf_genes <- c("IL6", "PDGFRA", "CXCL12", "CFD", "LMNA",
"AGTR1", "HAS1", "CXCL1", "CXCL2", "CCL2", "IL8")
mycaf_genes <- c("ACTA2", "TAGLN", "MMP11", "MYL9",
"HOPX", "POSTN", "TPM1", "TPM2")
pan_caf_genes <- c("COL1A1", "FAP", "PDPN", "DCN", "VIM")
gene_sets <- list(icaf_genes, mycaf_genes, pan_caf_genes)
names(gene_sets) <- c("iCAF", "myCAF", "Pan-CAF")
for (i in seq(gene_sets)) {
gene_sets[[i]] <- gene_sets[[i]][gene_sets[[i]] %in% rownames(fibro)]
}
fibro <- vamForSeurat(fibro,
gene.set.collection = gene_sets,
gamma = TRUE)
DefaultAssay(fibro) <- "VAMcdf"We can easily define cluster 0 as the myCAF population, and cluster 2 as the slightly smaller iCAF population. Cluster 4 shows no enrichment whatsoever for either the iCAF or myCAF gene sets.
p20 <- FeaturePlot(fibro, reduction = "tsne", features = "iCAF", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "iCAF") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p21 <- FeaturePlot(fibro, reduction = "tsne", features = "myCAF", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "myCAF") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p20 | p21) / p17So we have a small cluster of 23 cells that does not appear to express any of the fibroblast, CAF, perivascular, or endothelial markers. After performing differential expression analysis, we find that the top 3 markers for cluster 4 are CLU, CD74, and CRYAB. Interestingly, CLU and CD74 were found to be differentially expressed in the apCAF population discovered in the KPC mouse models of CAFs in Elyada et al.
DefaultAssay(fibro) <- "SCT"
fibro_markers <- FindAllMarkers(fibro,
assay = "SCT",
logfc.threshold = 1.5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE)
fibro_markers %>%
filter(cluster == 4) %>%
arrange(desc(avg_log2FC)) %>%
dplyr::select(cluster, gene, avg_log2FC, p_val_adj, pct.1, pct.2) %>%
slice_head(n = 5) %>%
kbl(booktabs = TRUE, digits = 4, row.names = FALSE) %>%
kable_minimal("hover", full_width = FALSE)| cluster | gene | avg_log2FC | p_val_adj | pct.1 | pct.2 |
|---|---|---|---|---|---|
| 4 | CLU | 4.1036 | 0 | 1.000 | 0.359 |
| 4 | CRYAB | 4.1022 | 0 | 0.957 | 0.365 |
| 4 | CD74 | 2.8741 | 0 | 0.957 | 0.317 |
| 4 | HLA-DRA | 2.5622 | 0 | 0.870 | 0.236 |
| 4 | HLA-DRB1 | 2.3438 | 0 | 0.609 | 0.149 |
We re-run GSEA, again using the VAM package and the differentially expressed genes for the apCAF population as defined in Elyada et al (with the mouse gene names converted to HGNC symbols). We can see that the apCAF pathway is strongly enriched in cluster 4 as compared to the other CAF clusters.
apcaf_genes <- c("HLA-DQB1", "CD74", "SAA3P", "SLPI")
gene_sets <- list(icaf_genes, mycaf_genes, apcaf_genes, pan_caf_genes)
names(gene_sets) <- c("iCAF", "myCAF", "apCAF", "Pan-CAF")
for (i in seq(gene_sets)) {
gene_sets[[i]] <- gene_sets[[i]][gene_sets[[i]] %in% rownames(fibro)]
}
fibro <- vamForSeurat(fibro,
gene.set.collection = gene_sets,
gamma = TRUE)
DefaultAssay(fibro) <- "VAMcdf"
p22 <- FeaturePlot(fibro, reduction = "tsne", features = "apCAF", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(title = "apCAF") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p22 / p17Finally, we add cell labels to our identified clusters.
fibro$label <- case_when(fibro$seurat_clusters == 0 ~ "myCAF",
fibro$seurat_clusters == 1 ~ "Perivascular",
fibro$seurat_clusters == 2 ~ "iCAF",
fibro$seurat_clusters == 3 ~ "Endothelial",
fibro$seurat_clusters == 4 ~ "apCAF")
Idents(fibro) <- "label"
p23a <- DimPlot(fibro, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p23aHere’s the marker genes for each cluster.
DefaultAssay(fibro) <- "SCT"
fibro_markers2 <- FindAllMarkers(fibro,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p23b <- DotPlot(fibro, features = fibro_markers2$gene, dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 12, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 16),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 18)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p23bGoing back to the main Seurat object, we should have a large population of T and NK cells that warrants further investigation. Using CD3D expression we can easily identify clusters 0, 3, and 7 as the mixed T & NK cells. We already see some good separation, so reclustering the cells should have good results.
p24 <- FeaturePlot(pdac, reduction = "tsne", features = "CD3D", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p24 / p2nkt_tumor <- subset(pdac, subset = seurat_clusters %in% c(0, 3, 8) & condition == "PDAC")Here we run ReclusterCells(), while treating the NK & T cell clusters as one large group. This will hopefully allow use to elucidate T cell subtypes. We use fewer PCs for these cells since the differences between immune cells are subtle, and adding more PCs will most likely only contribute noise to our analyses.
nkt_tumor <- ReclusterCells(nkt_tumor,
which.clust = list(0, 3, 8),
merge.clusters = TRUE,
n.HVG = 4000,
n.PC = 15,
k.vals = c(30, 40, 50, 60),
resolution.vals = c(.1, .2, .3),
redo.embedding = TRUE,
random.seed = 629)
nkt_tumor_pc <- Embeddings(nkt_tumor, "pca")## [1] "Reclustering cells in clusters 0, 3, 8 using k = 30 & resolution = 0.3; S = 0.455"
Once again we’ll run Fit-SNE on the reclustered cells.
# import data
nkt_tumor_pc = r.nkt_tumor_pc
# run Fit-SNE
affin_nkt_tumor = PerplexityBasedNN(nkt_tumor_pc, perplexity=100, metric='cosine', random_state=629)
init = initialization.pca(nkt_tumor_pc, random_state=629)
tsne_nkt_tumor = TSNEEmbedding(init, affin_nkt_tumor, negative_gradient_method='fft')
embed_nkt_tumor1 = tsne_nkt_tumor.optimize(n_iter=250, exaggeration=4, momentum=0.6)
embed_nkt_tumor2 = embed_nkt_tumor1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
affin_nkt_tumor.set_perplexity(50)
embed_nkt_tumor3 = embed_nkt_tumor2.optimize(n_iter=500, exaggeration=1, momentum=0.6)The reembedding isn’t perfect, which we somewhat expected as immune cells are difficult to tell apart based on the transcriptome alone.
embed_nkt_tumor <- as.matrix(py$embed_nkt_tumor3)
rownames(embed_nkt_tumor) <- colnames(nkt_tumor)
nkt_tumor@reductions$bh_tsne <- nkt_tumor@reductions$tsne
nkt_tumor@reductions$tsne<- CreateDimReducObject(embeddings = embed_nkt_tumor,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p25 <- DimPlot(nkt_tumor, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::category20_d3")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p25Clusters 0 and 1 contain our CD4+ T cells, which we characterize using IL7R as we did in the PBMC3k vignette. It’s somewhat outside of our scope here to determine which subtype each cluster belongs to, so we’ll simply assign both clusters the same label and move on.
p26 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "IL7R", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p27 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD69", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p26 | p27) / p25Cluster 3 contains the regulatory T cells.
p28 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "IL2RA", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p29 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "FOXP3", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p28 | p29) / p25We can find the proliferating T-regs in cluster 7.
p30 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "TOP2A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p30 / p25Mast cells can be identified using TPSAB1 expression in cluster 6.
p31 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "TPSAB1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p31 / p25NKG7 and PRF1 show us the NK cells in cluster 5.
p32 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "NKG7", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p33 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "PRF1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p32 | p33) / p25We use CD8A and CD2 to reveal the CD8+ T cells within clusters 2 and 4.
p34 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD8A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p35 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD2", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p34 | p35) / p25Lastly, we have cluster 8, which doesn’t highly express our pan-T/NK cell markers CD3D and CD2. It could be a myeloid cell cluster that was mistakenly grouped with the T/NK cells. We’ll start with a Wilcoxon test to determine its differentially expressed genes. Interestingly, several intermediate monocyte markers - LYZ, HLA-DRA, CD74, and HLA-DPB1 - are differentially expressed in cluster 8.
clust8_markers <- FindAllMarkers(nkt_tumor,
test.use = "wilcox",
min.diff.pct = .2,
logfc.threshold = .5,
verbose = FALSE,
random.seed = 629) %>%
filter(cluster == 8, p_val_adj < .05) %>%
arrange(desc(1 - p_val_adj))
clust8_markers %>%
filter(gene %in% c("LYZ", "HLA-DRA", "CD74", "HLA-DPB1")) %>%
dplyr::select(cluster, gene, avg_log2FC, p_val_adj, pct.1, pct.2) %>%
kbl(booktabs = TRUE, digits = 4, row.names = FALSE) %>%
kable_minimal("hover", full_width = FALSE)| cluster | gene | avg_log2FC | p_val_adj | pct.1 | pct.2 |
|---|---|---|---|---|---|
| 8 | HLA-DRA | 3.1725 | 0 | 0.974 | 0.410 |
| 8 | HLA-DPB1 | 2.2895 | 0 | 0.895 | 0.387 |
| 8 | CD74 | 2.7595 | 0 | 0.974 | 0.723 |
| 8 | LYZ | 1.9255 | 0 | 0.526 | 0.107 |
We’ll plot some of those markers below.
p36 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "LYZ", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p37 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "HLA-DRA", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p38 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "CD74", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p39 <- FeaturePlot(nkt_tumor, reduction = "tsne", features = "HLA-DPB1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
((p36 | p37) / (p38 | p39)) / p25We add labels to our T cell Seurat object and visualize the results.
nkt_tumor$label <- case_when(nkt_tumor$seurat_clusters == 0 ~ "CD4+ T",
nkt_tumor$seurat_clusters == 1 ~ "CD4+ T",
nkt_tumor$seurat_clusters == 2 ~ "CD8+ T",
nkt_tumor$seurat_clusters == 3 ~ "T-reg",
nkt_tumor$seurat_clusters == 4 ~ "CD8+ T",
nkt_tumor$seurat_clusters == 5 ~ "NK",
nkt_tumor$seurat_clusters == 6 ~ "Mast",
nkt_tumor$seurat_clusters == 7 ~ "Proliferating T-reg",
nkt_tumor$seurat_clusters == 8 ~ "Intermediate Monocyte")
Idents(nkt_tumor) <- "label"
p40a <- DimPlot(nkt_tumor, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 4)))
p40aHere’s the marker genes.
nkt_tumor_markers2 <- FindAllMarkers(nkt_tumor,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p40b <- DotPlot(nkt_tumor, features = unique(nkt_tumor_markers2$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 9, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 14)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p40bnkt_norm <- subset(pdac, subset = seurat_clusters %in% c(0, 3, 8) & condition == "AdjNorm")nkt_norm <- ReclusterCells(nkt_norm,
which.clust = list(0, 3, 8),
merge.clusters = TRUE,
n.HVG = 4000,
n.PC = 10,
k.vals = c(15, 20, 25),
resolution.vals = c(.2, .3, .4),
nn.metric = "euclidean",
redo.embedding = TRUE,
random.seed = 629)
nkt_norm_pc <- Embeddings(nkt_norm, "pca")## [1] "Reclustering cells in clusters 0, 3, 8 using k = 15 & resolution = 0.3; S = 0.438"
As with the other reclusterings, we’ll run Fit-SNE in order to (hopefully) obtain a better low-dimensional embedding of our cells.
# import data
nkt_norm_pc = r.nkt_norm_pc
# run Fit-SNE
affin_nkt_norm = PerplexityBasedNN(nkt_norm_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(nkt_norm_pc, random_state=629)
tsne_nkt_norm = TSNEEmbedding(init, affin_nkt_norm, negative_gradient_method='fft')
embed_nkt_norm1 = tsne_nkt_norm.optimize(n_iter=250, exaggeration=12, momentum=0.6)
embed_nkt_norm2 = embed_nkt_norm1.optimize(n_iter=750, exaggeration=1, momentum=0.8)embed_nkt_norm <- as.matrix(py$embed_nkt_norm2)
rownames(embed_nkt_norm) <- colnames(nkt_norm)
nkt_norm@reductions$bh_tsne <- nkt_norm@reductions$tsne
nkt_norm@reductions$tsne<- CreateDimReducObject(embeddings = embed_nkt_norm,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p41 <- DimPlot(nkt_norm, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p41First we’ll use a Wilcoxon test to determine which genes characterize each cluster.
nkt_norm_markers <- FindAllMarkers(nkt_norm,
logfc.threshold = .5,
min.diff.pct = .2,
verbose = FALSE,
only.pos = TRUE,
random.seed = 629) %>%
filter(p_val_adj < .05)
nkt_norm_markers %>%
dplyr::select(cluster, gene, avg_log2FC, p_val_adj, pct.1, pct.2) %>%
group_by(cluster) %>%
top_n(n = 3, wt = avg_log2FC) %>%
kbl(booktabs = TRUE, digits = 4) %>%
kable_minimal("hover", full_width = FALSE)| cluster | gene | avg_log2FC | p_val_adj | pct.1 | pct.2 |
|---|---|---|---|---|---|
| 0 | KLRB1 | 1.7367 | 0 | 0.843 | 0.316 |
| 0 | IL7R | 0.8630 | 0 | 0.897 | 0.612 |
| 1 | GZMK | 1.5317 | 0 | 0.925 | 0.350 |
| 1 | CST7 | 0.9069 | 0 | 0.846 | 0.452 |
| 1 | CCL4 | 1.0830 | 0 | 0.817 | 0.377 |
| 3 | KLRC1 | 0.8868 | 0 | 0.429 | 0.050 |
| 3 | GNLY | 1.2618 | 0 | 0.521 | 0.122 |
| 3 | HOPX | 0.6239 | 0 | 0.663 | 0.297 |
| 4 | GZMB | 2.3387 | 0 | 0.730 | 0.090 |
| 4 | NKG7 | 2.2958 | 0 | 0.873 | 0.515 |
| 4 | GNLY | 2.5414 | 0 | 0.611 | 0.149 |
| 5 | TNFRSF4 | 1.7814 | 0 | 0.691 | 0.043 |
| 5 | TNFRSF18 | 1.6798 | 0 | 0.660 | 0.047 |
| 5 | LTB | 1.7286 | 0 | 0.915 | 0.472 |
| 6 | SELL | 1.1920 | 0 | 0.761 | 0.080 |
| 6 | CCR7 | 1.0374 | 0 | 0.705 | 0.107 |
| 6 | EIF3E | 0.7694 | 0 | 0.830 | 0.520 |
| 7 | STMN1 | 2.7810 | 0 | 0.909 | 0.094 |
| 7 | TUBA1B | 3.1534 | 0 | 1.000 | 0.310 |
| 7 | HMGB2 | 3.1877 | 0 | 1.000 | 0.424 |
We can use IL7R and S100A4 expression to identify the memory CD4+ T cells in clusters 0 & 2. IL7R and CCR7 identify the naive CD4+ T cells in the adjacent cluster 6.
p42 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "IL7R", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p43 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "S100A4", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p44 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "CCR7", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p42 | p43 | p44) / p41Next we reveal the CD8+ T cells in clusters 1 and 3 with CD8A, as per usual.
p45 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "CD8A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p45 / p41The T-reg cluster, cluster 5, is identified using TIGIT and FOXP3.
p46 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "TIGIT", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p47 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "FOXP3", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p46 | p47) / p41The natural killers can be found in cluster 4 through their expression of PRF1 and NKG7. The NKG7 expression also confirms the identities of the CD8+ T cells we just annotated.
p48 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "PRF1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p49 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "NKG7", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p48 | p49) / p41The tiny proliferating T-reg population in cluster 7 is characterized by TOP2A.
p50 <- FeaturePlot(nkt_norm, reduction = "tsne", features = "TOP2A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p50 / p41We add final cell labels to our 8 cell clusters.
nkt_norm$label <- case_when(nkt_norm$seurat_clusters == 0 ~ "Memory CD4+ T",
nkt_norm$seurat_clusters == 1 ~ "CD8+ T",
nkt_norm$seurat_clusters == 2 ~ "Memory CD4+ T",
nkt_norm$seurat_clusters == 3 ~ "CD8+ T",
nkt_norm$seurat_clusters == 4 ~ "NK",
nkt_norm$seurat_clusters == 5 ~ "T-reg",
nkt_norm$seurat_clusters == 6 ~ "Naive CD4+ T",
nkt_norm$seurat_clusters == 7 ~ "Proliferating T-reg")
Idents(nkt_norm) <- "label"
p51a <- DimPlot(nkt_norm, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 4)))
p51aHere’s the marker genes.
nkt_norm_markers2 <- FindAllMarkers(nkt_norm,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p51b <- DotPlot(nkt_norm, features = unique(nkt_norm_markers2$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 10, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 14)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p51bWe use KRT8 expression to show the ductal cells residing in clusters 5 and 9. We note that some regions of clusters 5 and 9 have no KRT8 expression, which likely means that they are composed of different cell types.
p52 <- FeaturePlot(pdac, reduction = "tsne", features = "KRT8", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p52 / p2We run SCISSORS, then use a Wilcoxon differential expression test to determine potential markers for each cluster.
ductal <- ReclusterCells(pdac,
which.clust = list(5, 9),
merge.clusters = TRUE,
n.HVG = 4000,
n.PC = 20,
k.vals = c(20, 30, 40, 50), nn.metric = "euclidean",
resolution.vals = c(.2, .3, .4),
redo.embedding = TRUE,
random.seed = 629)
ductal_markers <- FindAllMarkers(ductal,
logfc.threshold = .5,
min.diff.pct = .2,
only.pos = TRUE,
verbose = FALSE,
random.seed = 629) %>%
filter(p_val_adj < .05)
ductal_markers %>%
dplyr::select(cluster, gene, avg_log2FC, p_val_adj, pct.1, pct.2) %>%
group_by(cluster) %>%
top_n(n = 3, wt = avg_log2FC) %>%
kbl(booktabs = TRUE, digits = 4) %>%
kable_minimal("hover", full_width = FALSE)## [1] "Reclustering cells in clusters 5, 9 using k = 50 & resolution = 0.2; S = 0.537"
| cluster | gene | avg_log2FC | p_val_adj | pct.1 | pct.2 |
|---|---|---|---|---|---|
| 0 | MUCL3 | 2.0586 | 0 | 0.823 | 0.204 |
| 0 | SLPI | 1.9466 | 0 | 0.976 | 0.698 |
| 0 | PSCA | 1.8493 | 0 | 0.705 | 0.208 |
| 1 | MUC1 | 1.0754 | 0 | 0.874 | 0.581 |
| 1 | PLCG2 | 1.5615 | 0 | 0.464 | 0.173 |
| 1 | LYZ | 1.5575 | 0 | 0.873 | 0.614 |
| 2 | CELA3A | 6.7301 | 0 | 0.796 | 0.053 |
| 2 | CTRB1 | 6.5305 | 0 | 0.741 | 0.058 |
| 2 | CTRB2 | 8.2024 | 0 | 0.847 | 0.125 |
| 3 | SERPING1 | 3.1176 | 0 | 0.996 | 0.266 |
| 3 | CLU | 3.5991 | 0 | 1.000 | 0.340 |
| 3 | SPP1 | 5.3722 | 0 | 0.992 | 0.416 |
| 4 | S100A2 | 4.4199 | 0 | 0.959 | 0.109 |
| 4 | CST1 | 3.6054 | 0 | 0.497 | 0.012 |
| 4 | PLAT | 3.2459 | 0 | 0.782 | 0.126 |
| 5 | APOA4 | 5.4298 | 0 | 0.698 | 0.002 |
| 5 | APOA1 | 5.8060 | 0 | 0.746 | 0.007 |
| 5 | ALDOB | 4.3565 | 0 | 0.968 | 0.043 |
| 6 | CRISP3 | 4.3529 | 0 | 0.930 | 0.022 |
| 6 | CRP | 3.6976 | 0 | 0.789 | 0.050 |
| 6 | TCN1 | 4.2364 | 0 | 0.895 | 0.119 |
Again, we’ll run Fit-SNE on the reclustered cells.
duct_pc <- Embeddings(ductal, reduction = "pca")# import data
duct_pc = r.duct_pc
# run Fit-SNE
affin_duct = PerplexityBasedNN(duct_pc, perplexity=30, random_state=629)
init = initialization.pca(duct_pc, random_state=629)
tsne_duct = TSNEEmbedding(init, affin_duct, negative_gradient_method='fft')
embed_d1 = tsne_duct.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_d2 = embed_d1.optimize(n_iter=750, exaggeration=1, momentum=0.8)We pull the results into R, making sure to save the Barnes-Hut t-SNE results in another reduction slot.
embed_duct <- as.matrix(py$embed_d2)
rownames(embed_duct) <- colnames(ductal)
ductal@reductions$bh_tsne <- ductal@reductions$tsne
ductal@reductions$tsne<- CreateDimReducObject(embeddings = embed_duct,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p53 <- DimPlot(ductal, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p53We use ANPEP & FABP1 expression to identify the lipid processing ductal cells in cluster 6.
p54 <- FeaturePlot(ductal, reduction = "tsne", features = "ANPEP", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p55 <- FeaturePlot(ductal, reduction = "tsne", features = "FABP1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p54 | p55) / p53Expression of SOD3 and CFTR reveals the secretory cells in cluster 3.
p56 <- FeaturePlot(ductal, reduction = "tsne", features = "SOD3", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p57 <- FeaturePlot(ductal, reduction = "tsne", features = "CFTR", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p56 | p57) / p53We use TFF1 and TFF2 expression to annotate the classical 1 epithelial cells in clusters 0 and 1. We can also see that the two classical 1 clusters are split by the sample from which the cells originate. Going forward, we’ll simply label both clusters as classical 1.
p58 <- FeaturePlot(ductal, reduction = "tsne", features = "TFF1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p59 <- FeaturePlot(ductal, reduction = "tsne", features = "TFF2", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p60 <- DimPlot(ductal, reduction = "tsne", group.by = "sample", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Sample") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p58 | p60 | p59) / p53The classical 2 cells are located in cluster 6, as evidenced by their expression of CRISP3.
p61 <- FeaturePlot(ductal, reduction = "tsne", features = "CRISP3", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p61 / p53The basal compartment weights from DECODER are highest in cluster 4, which we will denote as being composed of basal-like PDAC.
p62 <- FeaturePlot(ductal, reduction = "tsne", features = "basal") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Basal") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p62 / p53Lastly, we show that cluster 2 is composed of acinar cells using CTRB2.
p63 <- FeaturePlot(ductal, reduction = "tsne", features = "CTRB2", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p63 / p53We add final cluster labels to our Seurat object and visualize the results.
ductal$label <- case_when(ductal$seurat_clusters == 0 ~ "Classical 1",
ductal$seurat_clusters == 1 ~ "Classical 1",
ductal$seurat_clusters == 2 ~ "Acinar",
ductal$seurat_clusters == 3 ~ "Secretory",
ductal$seurat_clusters == 4 ~ "Basal",
ductal$seurat_clusters == 5 ~ "Lipid Proc.",
ductal$seurat_clusters == 6 ~ "Classical 2")
Idents(ductal) <- "label"
p64a <- DimPlot(ductal, reduction = "tsne",pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p64aHere’s the marker genes.
ductal_markers2 <- FindAllMarkers(ductal,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p64b <- DotPlot(ductal, features = unique(ductal_markers2$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 12, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 14)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p64bWe use JCHAIN (denoted IGJ in Elyada et al) to reveal the Plasma cells in cluster 10, and IRF7 to identify the plasmacytoid DCs in cluster 11.
p65 <- FeaturePlot(pdac, reduction = "tsne", features = "JCHAIN", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "JCHAIN (IGJ)") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p66 <- FeaturePlot(pdac, reduction = "tsne", features = "IRF7", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "IRF7") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p65 | p66) / p2We can identify the B cells in cluster 4 using joint expression of MS4A1 and CD79A. The cluster is split into two subclusters by tissue type: adjacent normal and PDAC. While it would be interesting to determine the genetic drivers of that separation, it’s somewhat outside of our scope here.
p67 <- FeaturePlot(pdac, reduction = "tsne", features = "MS4A1", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "MS4A1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p68 <- FeaturePlot(pdac, reduction = "tsne", features = "CD79A", pt.size = 0.75) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD79A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p69 <- DimPlot(pdac, reduction = "tsne", group.by = "condition", pt.size = 0.75) +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Tissue Type") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p67 | p69 | p68) / p2Lastly, we’ll split up the myeloid population by tissue condition (PDAC vs. adjacent normal) just like we did with the NK / T cells. We’ll run SCISSORS, annotate the clusters, and visualize the results.
First we’ll attempt to assign broad cell type labels to each of the four putative myeloid clusters. We’ll use the marker genes from Elyada et al once again.
myo_tumor <- subset(pdac, subset = seurat_clusters %in% c(1, 2, 7) & condition == "PDAC")
p70 <- DimPlot(myo_tumor, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p70Cluster 1 appears to be composed of our resident & alternatively activated macrophages due to its expression of CD14 & C1QA and SPP1, respectively.
p71 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "CD14", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD14") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p72 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "C1QA", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "C1QA") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p73 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "SPP1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "SPP1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p71 | p72 | p73) / p70We can use LYZ and S100A8 expression to reveal the classic monocytes and neutrophils in cluster 2.
p74 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "LYZ", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "LYZ") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p75 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "S100A8", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A8") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p74 | p75) / p70Lastly, we show that cluster 6 contains our various DC subtypes through its expression of FCER1A, a canonical dendritic cell marker.
p76 <- FeaturePlot(myo_tumor, reduction = "tsne", features = "FCER1A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "FCER1A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p76 / p70We’ll start with the Macrophages & Monocytes, since the DC populations are small and are best dealt with on their own.
myo_reclust <- ReclusterCells(myo_tumor,
which.clust = c(1, 2),
n.PC = 15,
merge.clusters = TRUE,
k.vals = c(40, 50, 60),
resolution.vals = c(.2, .3, .4),
n.HVG = 4000,
redo.embedding = TRUE,
random.seed = 629)
dc_reclust <- ReclusterCells(myo_tumor,
which.clust = 7,
n.PC = 15,
k.vals = c(20, 30, 40, 50),
resolution.vals = c(.3, .4, .5),
n.HVG = 4000,
nn.metric = "euclidean",
redo.embedding = TRUE,
random.seed = 629)## [1] "Reclustering cells in clusters 1, 2 using k = 60 & resolution = 0.2; S = 0.337"
## [1] "Reclustering cells in cluster 7 using k = 30 & resolution = 0.4; S = 0.429"
We’ll again run Fit-SNE on our reclustered cells.
mono_tumor_pc <- Embeddings(myo_reclust, "pca")
dc_tumor_pc <- Embeddings(dc_reclust, "pca")# import data
mono_pc = r.mono_tumor_pc
dc_pc = r.dc_tumor_pc
# Fit-SNE - monocytes & macrophages
affin_mono = PerplexityBasedNN(mono_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(mono_pc, random_state=629)
tsne_mono = TSNEEmbedding(init, affin_mono, negative_gradient_method='fft')
embed_mono1 = tsne_mono.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_mono2 = embed_mono1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
# Fit-SNE - DC
affin_dc = PerplexityBasedNN(dc_pc, perplexity=30, random_state=629)
init = initialization.pca(dc_pc, random_state=629)
tsne_dc = TSNEEmbedding(init, affin_dc, negative_gradient_method='fft')
embed_dc1 = tsne_dc.optimize(n_iter=250, exaggeration=12, momentum=0.6)
embed_dc2 = embed_dc1.optimize(n_iter=750, exaggeration=1, momentum=0.8)We pull the results back into R and visualize them.
embed_mono <- as.matrix(py$embed_mono2)
rownames(embed_mono) <- colnames(myo_reclust)
myo_reclust@reductions$bh_tsne <- myo_reclust@reductions$tsne
myo_reclust@reductions$tsne<- CreateDimReducObject(embeddings = embed_mono,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p77 <- DimPlot(myo_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p77embed_dc <- as.matrix(py$embed_dc2)
rownames(embed_dc) <- colnames(dc_reclust)
dc_reclust@reductions$bh_tsne <- dc_reclust@reductions$tsne
dc_reclust@reductions$tsne<- CreateDimReducObject(embeddings = embed_dc,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p78 <- DimPlot(dc_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p78First we ID the neutrophils in cluster 2 using S100A8 and S100A9 - marker genes used by Elyada et al.
p79 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "S100A8", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A8") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p80 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "S100A9", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A9") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p79 | p80) / p77We can identify the classical monocytes in cluster 1 through their expression of CD14 and lack of expression of CD16 aka FCGR3A, expression of which, alongside that of MS4A7, reveals the group of CD16+ monocytes in cluster 4. Finally, expression of those genes as well as S100A10 allows us to defined cluster 5 as containing intermediate monocytes.
p81 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "CD14", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD14") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p82 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "FCGR3A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "FCGR3A (CD16)") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p83 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "MS4A7", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "MS4A7") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p84 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "S100A10", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "S100A10") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p81 | p82 | p83 | p84) / p77Expression of C1QA, APOE, and SPP1 show us the resident and alternatively activated macrophages in clusters 0 and 3, respectively.
p85 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "C1QA", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "C1QA") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p86 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "APOE", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "APOE") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p87 <- FeaturePlot(myo_reclust, reduction = "tsne", features = "SPP1", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "SPP1") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p85 | p86 | p87) / p77We add cell labels to our Seurat object, then we’re off to the DCs.
myo_reclust$label <- case_when(myo_reclust$seurat_clusters == 0 ~ "Resident Macrophage",
myo_reclust$seurat_clusters == 1 ~ "Classical Monocyte",
myo_reclust$seurat_clusters == 2 ~ "Neutrophil",
myo_reclust$seurat_clusters == 3 ~ "Alt. Activated Macrophage",
myo_reclust$seurat_clusters == 4 ~ "CD16+ Monocyte",
myo_reclust$seurat_clusters == 5 ~ "Intermediate Monocyte")
Idents(myo_reclust) <- "label"We can use CLEC9A to annotate the cDC1 population in cluster 5.
p88 <- FeaturePlot(dc_reclust, reduction = "tsne", features = "CLEC9A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CLEC9A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p88 / p78High and low expression of CD1A and CD207 reveal the Langerhans-like DCB and DCA cells in clusters 3 and 4, respectively.
p89 <- FeaturePlot(dc_reclust, reduction = "tsne", features = "CD1A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD1A") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p90 <- FeaturePlot(dc_reclust, reduction = "tsne", features = "CD207", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "CD207") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p89 | p90) / p78Next we use LAMP3 to identify the activated DCs in cluster 6.
p91 <- FeaturePlot(dc_reclust, reduction = "tsne", features = "LAMP3", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "LAMP3") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p91 / p78Lastly, we use expression of two canonical cDC2 marker genes / transcription factors to identify clusters 0, 1, and 2 as cDC2 cells, which are split by sample ID.
p92 <- FeaturePlot(dc_reclust, reduction = "tsne", features = "CLEC10A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p93 <- FeaturePlot(dc_reclust, reduction = "tsne", features = "KLF4", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p94 <- DimPlot(dc_reclust, reduction = "tsne", group.by = "sample", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Sample") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p92 | p94 | p93) / p78We add final subcluster cell type labels to our Seurat object, then we’re off to the adjacent normal myeloid cells.
dc_reclust$label <- case_when(dc_reclust$seurat_clusters == 0 ~ "cDC2",
dc_reclust$seurat_clusters == 1 ~ "cDC2",
dc_reclust$seurat_clusters == 2 ~ "cDC2",
dc_reclust$seurat_clusters == 3 ~ "Langerhans-like DCB",
dc_reclust$seurat_clusters == 4 ~ "Langerhans-like DCA",
dc_reclust$seurat_clusters == 5 ~ "cDC1",
dc_reclust$seurat_clusters == 6 ~ "Activated DC")
Idents(dc_reclust) <- "label"p95a <- DimPlot(myo_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 4)))
p95aHere’s the marker genes.
myo_reclust_markers2 <- FindAllMarkers(myo_reclust,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p95b <- DotPlot(myo_reclust, features = unique(myo_reclust_markers2$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 9, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 14)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p95bp96a <- DimPlot(dc_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 2, override.aes = list(size = 3)))
p96aHere’s the marker genes.
dc_reclust_markers2 <- FindAllMarkers(dc_reclust,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p96b <- DotPlot(dc_reclust, features = unique(dc_reclust_markers2$gene), dot.scale = 8) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 9, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 14)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p96bWe have the same four clusters as in the tumor tissue - 1, 2, & 7 - and we’ll run the same analysis steps.
myo_norm <- subset(pdac, subset = seurat_clusters %in% c(1, 2, 7) & condition == "AdjNorm")myo_norm_reclust <- ReclusterCells(myo_norm,
which.clust = c(1, 2),
n.PC = 15,
merge.clusters = TRUE,
k.vals = c(20, 30, 40),
resolution.vals = c(.1, .2),
n.HVG = 4000,
redo.embedding = TRUE,
random.seed = 629)
dc_norm_reclust <- ReclusterCells(myo_norm,
which.clust = 7,
n.PC = 15,
k.vals = c(20, 30, 40, 50),
resolution.vals = c(.3, .4, .5),
n.HVG = 4000,
nn.metric = "euclidean",
redo.embedding = TRUE,
random.seed = 629)## [1] "Reclustering cells in clusters 1, 2 using k = 20 & resolution = 0.2; S = 0.328"
## [1] "Reclustering cells in cluster 7 using k = 30 & resolution = 0.3; S = 0.445"
For the last time, we run Fit-SNE in order to obtain a better embedding.
mono_norm_pc <- Embeddings(myo_norm_reclust, "pca")
dc_norm_pc <- Embeddings(dc_norm_reclust, "pca")# import data
mono_pc = r.mono_norm_pc
dc_pc = r.dc_norm_pc
# Fit-SNE - monocytes
affin_mono = PerplexityBasedNN(mono_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(mono_pc, random_state=629)
tsne_mono = TSNEEmbedding(init, affin_mono, negative_gradient_method='fft')
embed_mono1 = tsne_mono.optimize(n_iter=250, exaggeration=10, momentum=0.6)
embed_mono2 = embed_mono1.optimize(n_iter=750, exaggeration=1, momentum=0.8)
# Fit-SNE - DC
affin_dc = PerplexityBasedNN(dc_pc, perplexity=30, metric='cosine', random_state=629)
init = initialization.pca(dc_pc, random_state=629)
tsne_dc = TSNEEmbedding(init, affin_dc, negative_gradient_method='fft')
embed_dc1 = tsne_dc.optimize(n_iter=250, exaggeration=8, momentum=0.6)
embed_dc2 = embed_dc1.optimize(n_iter=750, exaggeration=1, momentum=0.8)We pull the results back into R and visualize them.
embed_mono <- as.matrix(py$embed_mono2)
rownames(embed_mono) <- colnames(myo_norm_reclust)
myo_norm_reclust@reductions$bh_tsne <- myo_norm_reclust@reductions$tsne
myo_norm_reclust@reductions$tsne<- CreateDimReducObject(embeddings = embed_mono,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p97 <- DimPlot(myo_norm_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p97embed_dc <- as.matrix(py$embed_dc2)
rownames(embed_dc) <- colnames(dc_norm_reclust)
dc_norm_reclust@reductions$bh_tsne <- dc_norm_reclust@reductions$tsne
dc_norm_reclust@reductions$tsne<- CreateDimReducObject(embeddings = embed_dc,
key = "FitSNE_",
assay = "SCT",
global = TRUE)
p98 <- DimPlot(dc_norm_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p98We use high S100A8 and S100A9 expression to define the neutrophils in cluster 4.
p99 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "S100A8", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p100 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "S100A9", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p99 | p100) / p97Next up are the classical monocytes in clusters 0 and 3, split by sample.
p101 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "LYZ", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p102 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "CD14", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p103 <- DimPlot(myo_norm_reclust, reduction = "tsne", group.by = "sample", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("miscpalettes::brightPastel")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2", title = "Sample") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p101 | p102 | p103) / p97C1QA and APOE show us the resident macrophages in clusters 1 and 2.
p104 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "C1QA", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p105 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "APOE", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p104 | p105) / p97Lastly, it seems we have a small group of CD8+ T cells in cluster 5 that snuck into the myeloid cluster, as defined by their expression of CD3D (marking them as NK / T cells), and CD8A & NKG7. I don’t believe they’re NK cells as they do not express PRF1 or GZMB.
p106 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "CD3D") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p107 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "CD8A") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p108 <- FeaturePlot(myo_norm_reclust, reduction = "tsne", features = "NKG7") +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p106 | p107 | p108) / p97We add labels to our clusters.
myo_norm_reclust$label <- case_when(myo_norm_reclust$seurat_clusters == 0 ~ "Classical Monocyte",
myo_norm_reclust$seurat_clusters == 1 ~ "Resident Macrophage",
myo_norm_reclust$seurat_clusters == 2 ~ "Resident Macrophage",
myo_norm_reclust$seurat_clusters == 3 ~ "Classical Monocyte",
myo_norm_reclust$seurat_clusters == 4 ~ "Neutrophil",
myo_norm_reclust$seurat_clusters == 5 ~ "CD8+ T")
Idents(myo_norm_reclust) <- "label"We annotate cluster 1 as conventional DC1 and cluster 0 as conventional DC2 through mutually exclusive expression of the canonical markers CLEC9A and CLEC10A, respectively.
p109 <- FeaturePlot(dc_norm_reclust, reduction = "tsne", features = "CLEC9A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
p110 <- FeaturePlot(dc_norm_reclust, reduction = "tsne", features = "CLEC10A", pt.size = 1.5) +
scale_color_gradientn(colors = wesanderson::wes_palette("Zissou1")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
NoLegend() +
theme(axis.title = element_blank())
(p109 | p110) / p98Labels are added to our adjacent normal tissue DC Seurat object.
dc_norm_reclust$label <- case_when(dc_norm_reclust$seurat_clusters == 0 ~ "cDC2",
dc_norm_reclust$seurat_clusters == 1 ~ "cDC2",
dc_norm_reclust$seurat_clusters == 2 ~ "cDC1")
Idents(dc_norm_reclust) <- "label"Here are the final annotations for the adjacent normal tissue myeloid population.
p111a <- DimPlot(myo_norm_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p111aAnd here’s their marker genes.
myo_norm_reclust_markers2 <- FindAllMarkers(myo_norm_reclust,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p111b <- DotPlot(myo_norm_reclust, features = unique(myo_norm_reclust_markers2$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 12, vjust = 0.5),
legend.position = "right", legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 14),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 14)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p111bHere’s the final DC annotations .
p112a <- DimPlot(dc_norm_reclust, reduction = "tsne", pt.size = 1.5) +
scale_color_manual(values = paletteer_d("ggsci::default_nejm")) +
labs(x = "Fit-SNE 1", y = "Fit-SNE 2") +
theme_yehlab() +
theme(plot.title = element_blank()) +
guides(color = guide_legend(nrow = 1, override.aes = list(size = 4)))
p112aAnd here are their marker genes.
dc_norm_reclust_markers2 <- FindAllMarkers(dc_norm_reclust,
logfc.threshold = .5,
test.use = "wilcox",
only.pos = TRUE,
random.seed = 629,
verbose = FALSE) %>%
filter(p_val_adj < .05) %>%
group_by(cluster) %>%
arrange(desc(avg_log2FC)) %>%
slice_head(n = 5)
p112b <- DotPlot(dc_norm_reclust, features = unique(dc_norm_reclust_markers2$gene), dot.scale = 10) +
scale_color_gradientn(colors = paletteer_d("wesanderson::Zissou1")) +
labs(color = "Expression", size = "% Expressed") +
theme(axis.text.x = element_text(angle = 90, size = 14, vjust = 0.5),
legend.position = "right",
legend.justification = "center",
panel.border = element_rect(fill = NA, size = 1, color = "black"),
axis.line = element_blank(),
legend.title = element_text(size = 16),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
axis.text.y = element_text(size = 16)) +
guides(color = guide_colorbar(title.position = "top",
barheight = unit(3, units = "cm"), title.hjust = 0.5),
size = guide_legend(title.position = "top", title.hjust = 0.5))
p112bWe’ll create a quick convenience function to help us save the figures.
SaveFigure <- function(my.plot = NULL, name = NULL, height = 8, width = 8) {
if (is.null(plot) | is.null(name)) stop("You forgot some arguments.")
# save figure as is - w/ axis labels, titles, etc.
dir <- "~/Desktop/R/SCISSORS/vignettes/figures_supp/Elyada"
ggsave(my.plot,
filename = paste0(name, ".pdf"),
device = "pdf",
units = "in",
path = dir,
height = height,
width = width)
# save "blank" figure w/ no labels, legends, etc.
dir <- "~/Desktop/R/SCISSORS/vignettes/figures_pub/Elyada"
plot_blank <- my.plot +
theme(axis.title = element_blank(),
panel.border = element_blank(),
plot.title = element_blank(),
plot.subtitle = element_blank(),
plot.caption = element_blank(),
legend.position = "none")
ggsave(plot_blank,
filename = paste0(name, ".pdf"),
device = "pdf",
units = "in",
path = dir,
height = height,
width = width)
}This section isn’t worth reading; it’s here solely to prove that the figures we present in our publication were dynamically generated during the knitting of this document.
SaveFigure(my.plot = p0a, name = "Seurat_Clusters_tSNE")
SaveFigure(my.plot = p0b, name = "Seurat_Clusters_Sil_Score")
SaveFigure(my.plot = p0c, name = "Seurat_Clusters_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p1, name = "Seurat_Clusters_UMAP")
SaveFigure(my.plot = p2, name = "Seurat_Clusters_FitSNE")
SaveFigure(my.plot = p3, name = "SingleR_Annos_SC_Ref")
SaveFigure(my.plot = p4, name = "SingleR_Annos_Bulk_Ref")
SaveFigure(my.plot = p5, name = "CONICSmat_Annos")
SaveFigure(my.plot = p6, name = "DECODER_Basal")
SaveFigure(my.plot = p7, name = "DECODER_Classical")
SaveFigure(my.plot = p8, name = "DECODER_Exocrine")
SaveFigure(my.plot = p9, name = "DECODER_Endocrine")
SaveFigure(my.plot = p10, name = "DECODER_Immune")
SaveFigure(my.plot = p11, name = "DECODER_Normal_Stroma")
SaveFigure(my.plot = p12, name = "DECODER_Act_Stroma")
SaveFigure(my.plot = p13, name = "Fibro_All_Cells_COL1A1")
SaveFigure(my.plot = p14, name = "Fibro_All_Cells_COL3A1")
SaveFigure(my.plot = p15, name = "Fibro_All_Cells_LUM")
SaveFigure(my.plot = p16, name = "Fibro_All_Cells_DCN")
SaveFigure(my.plot = p17, name = "Fibro_SCISSORS_FitSNE")
SaveFigure(my.plot = p18, name = "Fibro_Endothelial_PLVAP")
SaveFigure(my.plot = p19, name = "Fibro_Perivascular_RGS5")
SaveFigure(my.plot = p20, name = "Fibro_VAM_iCAF")
SaveFigure(my.plot = p21, name = "Fibro_VAM_myCAF")
SaveFigure(my.plot = p22, name = "Fibro_VAM_apCAF")
SaveFigure(my.plot = p23a, name = "Fibro_Final_Labels")
SaveFigure(my.plot = p23b, name = "Fibro_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p24, name = "NKT_All_Cells_CD3D")
SaveFigure(my.plot = p25, name = "NKT_Tumor_SCISSORS_FitSNE")
SaveFigure(my.plot = p26, name = "NKT_Tumor_CD4T_IL7R")
SaveFigure(my.plot = p27, name = "NKT_Tumor_CD4T_CD69")
SaveFigure(my.plot = p28, name = "NKT_Tumor_Treg_IL2RA")
SaveFigure(my.plot = p29, name = "NKT_Tumor_Treg_FOXP3")
SaveFigure(my.plot = p30, name = "NKT_Tumor_Prolif_Treg_TOP2A")
SaveFigure(my.plot = p31, name = "NKT_Tumor_Mast_TPSAB1")
SaveFigure(my.plot = p32, name = "NKT_Tumor_NK_NKG7")
SaveFigure(my.plot = p33, name = "NKT_Tumor_NK_PRF1")
SaveFigure(my.plot = p34, name = "NKT_Tumor_CD8T_CD8A")
SaveFigure(my.plot = p35, name = "NKT_Tumor_CD8T_CD2")
SaveFigure(my.plot = p36, name = "NKT_Tumor_Intermediate_Mono_LYZ")
SaveFigure(my.plot = p37, name = "NKT_Tumor_Intermediate_Mono_HLADRA")
SaveFigure(my.plot = p38, name = "NKT_Tumor_Intermediate_Mono_CD74")
SaveFigure(my.plot = p39, name = "NKT_Tumor_Intermediate_Mono_HLADPB1")
SaveFigure(my.plot = p40a, name = "NKT_Tumor_Final_Labels")
SaveFigure(my.plot = p40b, name = "NKT_Tumor_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p41, name = "NKT_Norm_SCISSORS_FitSNE")
SaveFigure(my.plot = p42, name = "NKT_Norm_CD4T_IL7R")
SaveFigure(my.plot = p43, name = "NKT_Norm_CD4T_S100A4")
SaveFigure(my.plot = p44, name = "NKT_Norm_CD4T_CCR7")
SaveFigure(my.plot = p45, name = "NKT_Norm_CD8T_CD8A")
SaveFigure(my.plot = p46, name = "NKT_Norm_Treg_TIGIT")
SaveFigure(my.plot = p47, name = "NKT_Norm_Treg_FOXP3")
SaveFigure(my.plot = p48, name = "NKT_Norm_NK_PRF1")
SaveFigure(my.plot = p49, name = "NKT_Norm_NK_NKG7")
SaveFigure(my.plot = p50, name = "NKT_Norm_Prolif_Treg_TOP2A")
SaveFigure(my.plot = p51a, name = "NKT_Norm_Final_Labels")
SaveFigure(my.plot = p51b, name = "NKT_Norm_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p52, name = "Ductal_All_Cells_KRT8")
SaveFigure(my.plot = p53, name = "Ductal_SCISSORS_FitSNE")
SaveFigure(my.plot = p54, name = "Ductal_Lipid_Proc_ANPEP")
SaveFigure(my.plot = p55, name = "Ductal_Lipid_Proc_FABP1")
SaveFigure(my.plot = p56, name = "Ductal_Secretory_SOD3")
SaveFigure(my.plot = p57, name = "Ductal_Secretory_CFTR")
SaveFigure(my.plot = p58, name = "Ductal_Classical1_TFF1")
SaveFigure(my.plot = p59, name = "Ductal_Classical1_TFF2")
SaveFigure(my.plot = p60, name = "Ductal_Classical1_SampleID")
SaveFigure(my.plot = p61, name = "Ductal_Classical2_CRISP3")
SaveFigure(my.plot = p62, name = "Ductal_DECODER_Basal")
SaveFigure(my.plot = p63, name = "Ductal_Acinar_CTRB2")
SaveFigure(my.plot = p64a, name = "Ductal_Final_Labels")
SaveFigure(my.plot = p64b, name = "Ductal_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p65, name = "Plasma_All_Cells_IGJ")
SaveFigure(my.plot = p66, name = "Plasmacytoid_DC_All_Cells_IRF7")
SaveFigure(my.plot = p67, name = "B_All_Cells_MS4A1")
SaveFigure(my.plot = p68, name = "B_All_Cells_CD79A")
SaveFigure(my.plot = p69, name = "B_All_Cells_Tissue_Type")
SaveFigure(my.plot = p70, name = "Myeloid_Tumor_All_Cells_FitSNE")
SaveFigure(my.plot = p71, name = "Myeloid_Tumor_All_Cells_CD14")
SaveFigure(my.plot = p72, name = "Myeloid_Tumor_All_Cells_C1QA")
SaveFigure(my.plot = p73, name = "Myeloid_Tumor_All_Cells_SPP1")
SaveFigure(my.plot = p74, name = "Myeloid_Tumor_All_Cells_LYZ")
SaveFigure(my.plot = p75, name = "Myeloid_Tumor_All_Cells_S100A8")
SaveFigure(my.plot = p76, name = "Myeloid_Tumor_All_Cells_FCER1A")
SaveFigure(my.plot = p77, name = "Myeloid_Tumor_SCISSORS_FitSNE")
SaveFigure(my.plot = p78, name = "DC_Tumor_SCISSORS_FitSNE")
SaveFigure(my.plot = p79, name = "Myeloid_Tumor_Neutrophil_S100A8")
SaveFigure(my.plot = p80, name = "Myeloid_Tumor_Neutrophil_S100A9")
SaveFigure(my.plot = p81, name = "Myeloid_Tumor_Classic_Mono_CD14")
SaveFigure(my.plot = p82, name = "Myeloid_Tumor_CD16_Mono_FCGR3A")
SaveFigure(my.plot = p83, name = "Myeloid_Tumor_CD16_Mono_MS4A7")
SaveFigure(my.plot = p84, name = "Myeloid_Tumor_CD16_Mono_S100A10")
SaveFigure(my.plot = p85, name = "Myeloid_Tumor_Resident_Macro_C1QA")
SaveFigure(my.plot = p86, name = "Myeloid_Tumor_Alt_Activated_Macro_APOE")
SaveFigure(my.plot = p87, name = "Myeloid_Tumor_Alt_Activated_Macro_SPP1")
SaveFigure(my.plot = p88, name = "DC_Tumor_cDC1_CLEC9A")
SaveFigure(my.plot = p89, name = "DC_Tumor_Langerhans_DC_CD1A")
SaveFigure(my.plot = p90, name = "DC_Tumor_Langerhans_DC_CD207")
SaveFigure(my.plot = p91, name = "DC_Tumor_Activated_DC_LAMP3")
SaveFigure(my.plot = p92, name = "DC_Tumor_cDC2_CLEC10A")
SaveFigure(my.plot = p93, name = "DC_Tumor_cDC2_KLF4")
SaveFigure(my.plot = p94, name = "DC_Tumor_Sample_ID")
SaveFigure(my.plot = p95a, name = "Myeloid_Tumor_Final_Labels")
SaveFigure(my.plot = p95b, name = "Myeloid_Tumor_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p96a, name = "DC_Tumor_Final_Labels")
SaveFigure(my.plot = p96b, name = "DC_Tumor_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p97, name = "Myeloid_Norm_SCISSORS_FitSNE")
SaveFigure(my.plot = p98, name = "DC_Norm_SCISSORS_FitSNE")
SaveFigure(my.plot = p99, name = "Myeloid_Norm_Neutrophil_S100A8")
SaveFigure(my.plot = p100, name = "Myeloid_Norm_Neutrophil_S100A9")
SaveFigure(my.plot = p101, name = "Myeloid_Norm_Classic_Mono_LYZ")
SaveFigure(my.plot = p102, name = "Myeloid_Norm_Classic_Mono_CD14")
SaveFigure(my.plot = p103, name = "Myeloid_Norm_Sample_ID")
SaveFigure(my.plot = p104, name = "Myeloid_Norm_Resident_Macro_C1QA")
SaveFigure(my.plot = p105, name = "Myeloid_Norm_Alt_Activated_Macro_APOE")
SaveFigure(my.plot = p106, name = "Myeloid_Norm_CD8T_CD3D")
SaveFigure(my.plot = p107, name = "Myeloid_Norm_CD8T_CD8A")
SaveFigure(my.plot = p108, name = "Myeloid_Norm_CD8T_NKG7")
SaveFigure(my.plot = p109, name = "DC_Norm_cDC1_CLEC9A")
SaveFigure(my.plot = p110, name = "DC_Norm_cDC2_CLEC10A")
SaveFigure(my.plot = p111a, name = "Myeloid_Norm_Final_Labels")
SaveFigure(my.plot = p111b, name = "Myeloid_Norm_Dotplot", height = 8, width = 12)
SaveFigure(my.plot = p112a, name = "DC_Norm_Final_Labels")
SaveFigure(my.plot = p112b, name = "DC_Norm_Dotplot", height = 8, width = 12)And of course:
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] celldex_1.0.0 wesanderson_0.3.6
## [3] kableExtra_1.3.4 reticulate_1.18
## [5] CONICSmat_0.0.0.1 paletteer_1.3.0
## [7] latex2exp_0.5.0 patchwork_1.1.1
## [9] mixtools_1.2.0 SCISSORS_0.0.2.0
## [11] SingleCellExperiment_1.12.0 data.table_1.14.0
## [13] cluster_2.1.1 biomaRt_2.46.3
## [15] decoderr_0.0.0.9000 SingleR_1.4.1
## [17] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [19] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7
## [21] IRanges_2.24.1 S4Vectors_0.28.1
## [23] BiocGenerics_0.36.0 MatrixGenerics_1.2.1
## [25] matrixStats_0.58.0 ggplot2_3.3.3
## [27] SeuratObject_4.0.0 Seurat_4.0.1
## [29] dplyr_1.0.5 VAM_0.5.2
## [31] Matrix_1.3-2 MASS_7.3-53.1
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.1 tidyselect_1.1.0
## [3] RSQLite_2.2.6 AnnotationDbi_1.52.0
## [5] htmlwidgets_1.5.3 grid_4.0.4
## [7] BiocParallel_1.24.1 Rtsne_0.15
## [9] munsell_0.5.0 codetools_0.2-18
## [11] ica_1.0-2 statmod_1.4.35
## [13] scran_1.18.6 future_1.21.0
## [15] miniUI_0.1.1.1 withr_2.4.2
## [17] colorspace_2.0-0 highr_0.8
## [19] knitr_1.32 rstudioapi_0.13
## [21] ROCR_1.0-11 tensor_1.5
## [23] listenv_0.8.0 labeling_0.4.2
## [25] GenomeInfoDbData_1.2.4 polyclip_1.10-0
## [27] pheatmap_1.0.12 farver_2.1.0
## [29] bit64_4.0.5 parallelly_1.24.0
## [31] vctrs_0.3.7 generics_0.1.0
## [33] xfun_0.22 BiocFileCache_1.14.0
## [35] squash_1.0.9 R6_2.5.0
## [37] phateR_1.0.7 rsvd_1.0.3
## [39] locfit_1.5-9.4 bitops_1.0-6
## [41] spatstat.utils_2.1-0 cachem_1.0.4
## [43] DelayedArray_0.16.3 assertthat_0.2.1
## [45] promises_1.2.0.1 scales_1.1.1
## [47] gtable_0.3.0 beachmat_2.6.4
## [49] globals_0.14.0 goftest_1.2-2
## [51] rlang_0.4.10 systemfonts_1.0.1
## [53] splines_4.0.4 lazyeval_0.2.2
## [55] prismatic_1.0.0 spatstat.geom_2.1-0
## [57] BiocManager_1.30.12 yaml_2.2.1
## [59] reshape2_1.4.4 abind_1.4-5
## [61] httpuv_1.5.5 tools_4.0.4
## [63] ellipsis_0.3.1 spatstat.core_2.0-0
## [65] jquerylib_0.1.3 RColorBrewer_1.1-2
## [67] ggridges_0.5.3 Rcpp_1.0.6
## [69] plyr_1.8.6 sparseMatrixStats_1.2.1
## [71] progress_1.2.2 zlibbioc_1.36.0
## [73] purrr_0.3.4 RCurl_1.98-1.3
## [75] prettyunits_1.1.1 rpart_4.1-15
## [77] openssl_1.4.3 deldir_0.2-10
## [79] pbapply_1.4-3 cowplot_1.1.1
## [81] zoo_1.8-9 ggrepel_0.9.1
## [83] magrittr_2.0.1 RSpectra_0.16-0
## [85] scattermore_0.7 lmtest_0.9-38
## [87] RANN_2.6.1 fitdistrplus_1.1-3
## [89] hms_1.0.0 mime_0.10
## [91] evaluate_0.14 xtable_1.8-4
## [93] XML_3.99-0.6 gridExtra_2.3
## [95] compiler_4.0.4 tibble_3.1.1
## [97] KernSmooth_2.23-18 crayon_1.4.1
## [99] htmltools_0.5.1.1 segmented_1.3-3
## [101] mgcv_1.8-34 later_1.1.0.1
## [103] tidyr_1.1.3 DBI_1.1.1
## [105] ExperimentHub_1.16.0 dbplyr_2.1.1
## [107] rappdirs_0.3.3 igraph_1.2.6
## [109] pkgconfig_2.0.3 scuttle_1.0.4
## [111] plotly_4.9.3 spatstat.sparse_2.0-0
## [113] xml2_1.3.2 svglite_2.0.0
## [115] bslib_0.2.4 dqrng_0.2.1
## [117] webshot_0.5.2 XVector_0.30.0
## [119] rvest_1.0.0 stringr_1.4.0
## [121] digest_0.6.27 sctransform_0.3.2
## [123] RcppAnnoy_0.0.18 spatstat.data_2.1-0
## [125] rmarkdown_2.7 leiden_0.3.7
## [127] edgeR_3.32.1 uwot_0.1.10
## [129] DelayedMatrixStats_1.12.3 curl_4.3
## [131] kernlab_0.9-29 shiny_1.6.0
## [133] lifecycle_1.0.0 nlme_3.1-152
## [135] jsonlite_1.7.2 BiocNeighbors_1.8.2
## [137] limma_3.46.0 viridisLite_0.4.0
## [139] askpass_1.1 fansi_0.4.2
## [141] pillar_1.6.0 lattice_0.20-41
## [143] fastmap_1.1.0 httr_1.4.2
## [145] survival_3.2-10 interactiveDisplayBase_1.28.0
## [147] glue_1.4.2 png_0.1-7
## [149] BiocVersion_3.12.0 bluster_1.0.0
## [151] bit_4.0.4 nnls_1.4
## [153] stringi_1.5.3 sass_0.3.1
## [155] rematch2_2.1.2 blob_1.2.1
## [157] AnnotationHub_2.22.0 BiocSingular_1.6.0
## [159] memoise_2.0.0 irlba_2.3.3
## [161] future.apply_1.7.0